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ABSTRACT 

In this paper, we present a new implementation of feedback due to active galactic nuclei 
(AGN) in cosmological simulations of galaxy formation. We assume that a fraction of jet 
energy, which is generated by an AGN, is transferred to the surrounding gas as thermal en- 
ergy. Combining a theoretical model of mass accretion onto black holes with a multiphase 
description of star-forming gas, we self-consistently follow evolution of both galaxies and 
their central black holes. The novelty in our model is that we consider two distinct accretion 
modes: standard radiatively efficient thin accretion disks and radiatively inefficient accretion 
flows which we will generically refer to as RIAFs; motivated by theoretical models for jet pro- 
duction in accretion disks, we assume that only the RIAF is responsible for the AGN feedback. 
The focus of this paper is to investigate the interplay between galaxies and their central black 
holes during the formation of a disc galaxy. We find that, after an initial episode of bursting 
star formation, the accretion rate onto the central black hole drops so that the accretion disk 
switches to a RIAF structure. At this point, the feedback from the AGN becomes efficient 
and slightly suppresses star formation in the galactic disk and almost completely halts star 
formation in the bulge. This suppression of the star formation regulates mass accretion onto 
the black hole and associated AGN feedback. As a result, the nucleus becomes a stochas- 
tically fuelled low-luminosity AGN (Seyfert galaxy) with recurrent short-lived episodes of 
activity after the star bursts. During the "on" events the AGN produces reasonably power- 
ful jets (radio-loud state) and is less luminous than the host galaxy, while in the "off" phase 
the nucleus is inactive and "radio-quiet". Our model predicts several properties of the low- 
luminosity AGN including the bolometric luminosity, jet powers, the effect on kpc-scale of 
the radio jet and the AGN lifetime, which are in broad agreement with observations of Seyfert 
galaxies and their radio activity. We also find that the ratios between the central black hole 
mass and the mass of the host spheroid at z = are ~ 10~ 3 regardless of the strength of 
either supernova feedback or AGN feedback because the radiation drag model directly relates 
the star formation activity in the galactic centre and the mass accretion rate onto the central 
black hole. 
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1 INTRODUCTION 

Recent observations have suggested a fundamental connec- 
tion between active galactic nuclei (AGN) and the formation 
of galaxies. Firstly, there is now a well established correla- 
tion between the properties of galaxies and the masses of 
the black holes (BH) at their centres. The mass of the cen- 
tral black holes tightly correlates with the mass of the host 
galactic bulges (with a median BH-to-bulg e mass ratio of 
~ .0 01, e.g. Kormendv & Richstond 1995 : Magorrian et al.l 
1 19981 : iMerritt & Ferraresd l200ll : iMcLure & Dunlod |2002| ; 
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iMarconi & Huntll2003l). as well as with the stellar velocity disper- 
sion of the bulges dFerrarese & MerritfcOOOl : iGebhardt et alB oOO ; 



iTremaine et alj l2002h . This discovery points to a fundamental 
connection between the growth of central BHs and the formation 
of stellar spheroids in galaxies. Secondly, high resolution X-ray 
observations of galaxy clusters have revealed large, radio-plasma 
cavities in intra-cluster medium. These are usually associated with 
episodic outbursts from a central radio galaxy and indicate that 
huge amounts of mechanical energy are being deposi ted into the 
intracluster medium by powerful AGN-driven jets (e.g.lAllen et al.l 
20061: iBirzan et alj|2004l : iFabian et alj|2006l: iRaffertv et ai]|2006l: 



Taylor et al.l I2006D . Simulations of the impact of these AGN- 



driven radio cavities suggest that this powerful feedback provides 
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sufficient energy to offset the cooling radiation from the cluster 
and potentially explain why so littl e cool (T < 3 keV) gas is 
seen in these systems (Ouilisetal. 2001; Churazov et al. 2002; 



Dalla Vecchia et al l 120041 : lOmma et all 120041 : ISiiacki & Springell 
2006l ; ISiiackietal.ll2007» . 



Motivated by these discoveries, simple prescriptions for the 
feedback from AGN (aka. "radio-mode feedback") have recently 
bee n incorporated into semi-analytic galaxy formation models 
(e.g. Granato et al. 2004; Monaco & Fontanot 2005; Cattaneo et al. 



l2005l ; ICroton et al .11 20061 ; foe Lucia et al]|2006l : lBower et alj|2006h . 
Including radio-mode AGN feedback has resulted in dramatic im- 
provements in the models' ability to match the sharp decline of the 
galaxy luminosity function and to explain the "down- sizing" seen 
in the evolution of the galaxy population. In particular. lBower et al.l 
shows that, by assuming that AGN feedback operates only in quasi- 
hydrosta tic halos where the cooling time is longer than the dynami- 
cal time dDalla Vecchia 2005; Sii acki & Spr ingel 2006), the galaxy 
luminosity fu nctions in local and hig her redshifts universe can be 
matched well. lCroton et alj d2006h and lKitzbichler & White! d2007h 
showed that similar results were obtained b y modifying the Bondi- 
Hoyle-Lyttleton accretion rate form ula dHovle & Lvtfletonlll939l ; 
Bondi & Hovle|[l944l ; [Bondi 195^) to account for the multiphase 
structure of the accreting gas in rapidly cooling halos. 

Recent simulations also have begun to track the impact of 
AGN feedback on t he gal a xy population dKawata & Gibsonll 20051; 



Siiacki & Spring ell 20061; |Di Matteo et al. 2005J; ISpringel et al.1 
2005bl ; ISiiacki et al.ll2007l : foi Matteo et al.ll2007l) . One of the first 
hydrodynamic simulation of g alaxy formation wh i ch in voked AGN 
feedback was performed by iKawata & Gibsorj (2005) by using 
cosmological simulations with smoothed particle hydrodynamics 
(SPH). They reproduced observed X-ray and optical properties of 
elliptical galaxies by injecting thermal energy into the centre of the 
main progenitor at z < 1, assuming that sufficient BH mass was 
present for the AGN to be active when a convergent gas inflow ex- 
ists at the centre of the main progenitor. A self-consistent model 
for AGN "quasar-mode" feedback associated wi th BH accretion 
in sim ulati ons of galaxy mergers w as proposed by fol Matte o et al.1 
d2005h and lSpringel et alj (2005b). They estimated accretion rate 
onto BHs by using a Bondi-Hoyle-Lyttleton parameterisation and 
injected some fraction of accreted rest mass energy into the sur- 
rounding interstellar medium (ISM) in the form of thermal energy. 
By using this model and performing a series of merger simulations, 
iHopkins et al. (2006) simulated evolution of quasar luminosity and 
predicted luminosity functions of quasars. 

In this paper, we introduce a new methodology to incorporate 
BH growth and AGN radio feedback self-consistently in cosmolog- 
ical simulations of galaxy formation. We are motivated to do this 
by the recent semi-analytic models that we have discussed above. 
These indicate that the distinction between "quasar" and "radio" 
modes is necessary to achieve good matches to the global proper- 
ties of galaxies. We distinguish these two modes depending on the 
accretion rate onto the central BHs assuming that the radio-mode 
exists only when the accretion rate onto the BH is low. 

We estimate the mass accretion rat e onto the central BHs 
by co mbining a radiation drag model by Kawakatu & Umemural 
d2002l) with our description of the ISM and star formation. A virtue 
of the radiation drag model is that it relates star formation activity 
in galactic centre, which can be resolved in our simulations, to mass 
accretion onto a central BHs. It implies that a fixed fraction of the 
star formation in bulges is convert ed into BH mass just as assumed 
in mo st semi-analytic models (e.g. lCroton et al.| [2006; Bow er et"al] 
2006). Since our estimate of the accretion rate onto central BHs is 



based on radiation from stars, the mass growth of BHs in our sim- 
ulations is resolution independent so long as the star formation is 
modelled in a resolution independent fashion. 

In order to model the AGN radio-mode feedback, we distin- 
guish two fundamentally different modes of AGN accretion: radia- 
tively efficient geometrica lly thin accretion flows (standard disks; 
Shakura & Sunvaev|[l973h and geometrically thick, radiatively in- 
efficient ^_accrrtion_Jtowj>_(whi£h_w^ 
RIAFs; 



accretion Hows (which we will genencall y reter to as 
Naravan et al. I I 19981 ; iNaravai] 12005; Ne mmen et al1l2006l; 



al. 
5). 



lYuanl2007l) . We assume only the latter are responsible fo r the radio- 
mode feedback through production of powerful jets (e.g.lRees et 
1 19821; lMeierll200ll ; iMaccarone et"ai1l2003l: IChurazov et al .1120051 
This assertion is supported by theoretical models for jet production 
since the jet power is intimately linked to the scale of the magnetic 
field loops in the vicinity of the last stable orbit aroun d the BH (e.g. 
iLivio et aljfl999l: lMeieij|200ll ; iNemmen et ai]|2007l) . Since RIAFs 
exist o nly when the accret ion rate is much lower than the Eddington 
rate dNaravan et aT1ll998l) , AGNs can produce strong feedback in 
objects which have low specific star formation rates because the ra- 
diation drag model directly connects the star formation rate around 
BHs and mass accretion rate onto the BHs. Given the limited nu- 
merical resolution inherent in our simulations, we adopt a subgrid 
model for the impact of the AGN feedback. We simply assume that 
a certain fraction of the available jet power thermally couples to 
nearby diffuse halo gas. 

Until today, most of simulations of galaxy formation with 
AGN feedback aimed to form elliptical galaxies because the AGN 
feedback is the most plausible mecha nism that halts the star 
formation activity in elliptical galaxie s (Kawata & Gibson 2005; 
ISpringel et"al1 l2005bl : foi Matteo etal] 120051 ; ISiiacki et al.N2007t) . 
Therefore roles of AGN feedback in disk galaxy formation are 
still largely unknown in spite that disk galaxies also harbour cen- 
tral supermassive BHs in their bulges. In this paper we hence fo- 
cus on the impact that feedback from an AGN jet might have on 
the formation of a large disk galaxy, and the resulting coevolution 
of the AGN and the host galaxy. The majority of AGN hosted by 
disk galaxies are known to be either Seyfert nuclei or the lowest 
luminosity AGN, the low- ioni sation nu clear emission-line regions 
(LINERs) dHo et al.|[l997l ; see lHoll2004l for a r eview). Several tar- 
geted observations of individual objects (e.g., IWilson & Ulvestaa 
1 19831; iMorganti. Oosterloo & Tsvetanovlll998l; ICecil et alJbOobT 
Whittle & Wilson! |2004| ; IKeel et all 120061 : iKharb et al.1 I2OO6I; 
Middelberg et al. 2007) indicate that Seyfert galaxies commonly 
host compact radio jets spanning ~ 10 — 100 pc which strongly 
affect the st ate of the ISM and contribute to its ionisation via 
shocks (e.gjKukula et all! 19931, 1 1 9961: iFalcke Wilson. & Simpsor] 



ll998l;ICapetti et alj 19991 ; iBicknell et alj 19981 ; iMundell et alj2003l; 
Riffel et al. 2006). It is also common that kpc-scale radio structures 
are o bserved, which are though t as extensions of the inner radio jets 
(e.g. IWilson & Ulvestadll 19831; IMorganti, Oosterloo. & Tsvetanoyl 
19981; ICecil et alj|2000l: IWhittle & Wilsorj|2004l: IKeel et al.ll2006l: 



Kharb et alj|2006l ; iMiddelberg et alj I2007T) , and so metimes blow 
out of the galactic disk with n o preferential direction dSchmitt et al.1 
120011 ; IGallimore et al.ll200o1) . Very Large Array surveys of sam- 
ples of Seyfert and LINER galaxies indicate that radio outflows 
with extent > 1 kpc are a common feature in Seyfert galax- 
ies, and presumably are dri v en by the AGN jets dColbert et al] 
1 19961 : IGallimore etal J l2006h . IGallimore et alj d2006l) in particu- 
lar find that > 44% of the galaxies in their complete sample 
display kpc-scale radio outflows, and they cannot rule out that 
most Seyfert galaxies produce extended radio outflows. The ki- 
netic power carried by the Seyfert jets as estimated from obser- 
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vations can be as high as 10 43 erg s 1 (e.g., Iwilson & Ulvestad 



1 19831 : iMorganti. Oosterloo. & Tsvetanov|[l998l : lKharb et alj|2006t) . 
If an appreciable fraction of this power reaches out of the nuclear 
region, then the jets may affect the environment of the AGN and be 
a source of feedback in Seyfert galaxies. The impact of feedback 
from AGN jets on the nuclear environment of the host disk galax- 
ies and the evolution of the nuclear activity during the formation of 
the disk galaxy are open problems, which we attack using our cos- 
mological simulations. In this paper, we show that distinguishing 
the two modes of AGN accretion yields results which are in broad 
agreement with observations of Seyferts. 

The paper is organised as follows. We first describe our micro- 
scopic descriptions of the ISM, star formation, and stellar feedback 
in Section 2 and BH growth and AGN feedback in Section 3. We 
present the details of our cosmological simulations in Section 4 and 
the results in Section 5. Finally, we summarise and discuss our re- 
sults in Section 6. Readers who are observation-oriented and are 
not interested in the technical details of the simulation may wish 
to skip directly to Section [6] where we compare our results with 
observations of nuclear activity in Seyfert galaxies. 



2 THE MODEL OF STAR-FORMING GAS AND 
STELLAR FEEDBACK 

In order to study the effect of AGN feedback on galaxy for- 
mation, we construct a physically motivated, well-controlled 
model of the star-forming gas. Based on the simplest star for- 
mation and feedback recipes used in the first generation of hy- 
drodynamic simulations of galaxy formation (N avarro & White! 



1993 JSteinmetz & Muellerll 994; Cen & Ostriker 1992; Katz et al. 



1996), several authors have introduced 'multiphase' models for 
star-forming gas in which the ISM is treated as a number of distinct 
phases by formulating differ ential equations th at model the inter- 
actions between the phases I Yepes et al] 1997 ; S pringe l & Hern- 
quist 2003 hereafter ISHOl lsamland & Gerhard 12003 ; IOEFJ05I; 
IScannapieco et al]|2006l but see iBooth et al.ll2007l where they ex- 
plicitly model 'clouds' by means of sticky parti cles in SP H sim- 
ulations). Our model is based on the one used in OEFJ05 but has 
been significantly modified. We describe the details in the follow- 
ing subsections. 



2.1 Gas cooling 

We compute the non-equilibrium cooling/heating rate and ionisa- 
tion state of each g as particle by using the cooling functions in 
iTheuns et af] J2OO2K which include cooling by H, He, and metals, 
in the presence of an e volving but uniform ultraviolet background 
( Haard t & Madav] 1 19961) that is switched on at 2 ~ 6. Inverse 
Compton cooling is also included. 

Metals are carried by particles and once assigned to a parti- 
cle, they remain with it. Nonetheless, effective mixing takes place 
because we use the smoothed metallicity (smoothed in the same 
way as other SPH quantities) when computing the cooling rate and 
defining the metallicity of a newly-formed star. 



2.2 Cloud formation and the cloud model 

As in ISH03l andlOEFJ05 we treat an SPH particle as a hybrid par- 
ticle that consists of two distinct phases, i.e. hot ambient gas and 
cold clouds, once its gas density, p, exceeds a density threshold, 
Pth- 



We assume that the cold clouds form and grow through ther- 
mal instability, that is 



dM c 



dM h 



dt 



A net (ph,iih, Z) M h 



Ph 



(1) 



TI 

where M c and Af n are mass in cold clouds and mass in the hot 
phase associated with the particle, respectively, A nct is the cooling 
function for gas of metallicity Z, and Uh and u c represent specific 
internal energy of the hot phase and cold clouds, respectively. This 
implies that we assume that the gas is thermally unstable when p > 
pth- Cold clouds remain at a fixed temperature, T c — WOK hence 
u c is a constant as well. Note, however, as far as T c is much lower 
than 10 4 K our results are hardly affected by the choice of T c . 
We allow cold clouds to have a mass spectrum 



, . _ dJVc 
(m) = — — oc m 
dm 



(2) 



Observationally a is in a range between a = 1.5 to 1. 9 
dSolomon & Rivolol Il989l ; iFukui et al.l l200ll; iHever et al.l l200ll) . 
Since the effect of changing the value of a is completely degen- 
erate with changing other parameters we will introduce such as star 
formation efficiency and thermal conduction efficiency, we simply 
fix a — 1.7 and the mass range of clouds, 1O 2 -1O 6 M0, throughout 

this paper. 

Following Sam land & Gerhard! §003), we assume spherical 
clouds which follow the mass-radius relation o f lElmegreerJd 19891) 



where P4 = 3 is the pressure in the ambient hot phase 

which is equal to the effective pressure of the particle. 



2.3 Star formation 

We assume that giant 

to 1O 6 M are eligible to form stars with a star formation timescale 
that is proportional to the dynamical time of each cloud. Therefore 
the star formation rate for a cloud of mass m is 
dm* 



m m 

dt f s f tdyn 



(4) 



where m, is stellar mass, t s t = fdyn/c* is the star formation 
timescale for a cloud of mass m, c* is a dimensionless star for- 
mation efficiency parameter that must be less than unity, and td yn 
is the dynamical timescale of the cloud. For the mass-radius rela- 
tion in equation y), the dynamical timescale of a cloud of mass m 
is given as 



tdyn 



3tt 



-1 1/2 



32Gp(m) 



0.32 P A 



-3/8 



m 



1/4 



Myr, 



(5) 



where p(m) is the mean density of a cloud of mass m. 

It should be noted that now the star formation timescale is a 
function of ambient pressure and it could be very short in high- 
pressure mediu m. The refore we naturally have 'shock-induced' 
starbursts that lOEFJ05l introduced by changing the star formation 
efficiency by hand when the rate of change of the entropy by shock 
heating exceeds a threshold value. The total star formation rate for 
clouds of total mass of 1 Mq with the mass spectrum (f>(m) is ob- 
tained by integrating m„ over the range of masses that are eligible 
to form stars, 



rW°M Q 

M„ = / m*(m)<^(m)dm, 

J 1() 4 M 



(6) 
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where the normalisation 

10 6 Af Q 
1O 2 M 

is applied. 



m<f)(m)dm — 1 



2.4 Cloud evaporation 

The differe nt phases of the ISM exch ange mass by evapora- 
tion. As in ISamland & Gerhard! d2003h . we use the model of 
ICowie & McKeej dl977h for the cloud evaporation, that is, 



7TIEVP 



dm 
~dt 



4.4 x 10" 



167r finr 
25fc B 



??EVP 



T 



i/2 



m 



Mq Myr 



-1/4^,5/2 



(7) 



where ?7evp, p, k, r, and ftp, are the dimensionless conduction ef- 
ficiency parameter, the mean m olecular weigh t, the conductivity 
where we use the Spitzer value dSpitzei|[T96a) . the radius of the 
cloud, and the Boltzmann constant. The total evaporation rate for 
clouds of total mass 1 Mq then becomes 



Me 



1O 6 M 
1O 2 M 



rhEVp(rn)(j>(m)drn. 



(8) 



2.5 Supernova feedback and chemical enrichment 

We consider each stellar particle as a single stellar population (SSP) 
having its own age, metallicity and initial mass function (IMF) 
(j>„(m). Through out this paper w e assume that the IMF is always 
the Salpeter IMF l Sal peteJl955h whose lower and upper mass lim- 
its are 0.1 and 100 M©, respectively. The recycled mass fraction, 
type II and type la supernova (SN) rates, and metal production rates 
for some species are calculated according to the age, metallicity, 
and IMF of the SSP. Each SN is assumed to give t^snIO 51 erg of 
energy to surrounding gas, where ?7sn is t he feedback efficiency 
parameter. We refer the readers to OEFJ05 and references therein 
for a more detailed description. 



2.6 Numerical implementation 

We have implemented the physics descri bed above into a publicly 
available PM-TreeSPH code GADGE T2 dSpringell2005h. a succes- 
sor of the TreeSPH code GADGET dSpringel et alj|200lh . We do 
not use the instantaneous recycling approximation (IRA) which is 
often assumed in studies of galaxy formation. In this section we de- 
scribe how the physical processes described above are implemented 
in our code. 

For p > pth, an SPH particle is eligible to form cold clouds. 
Hence it is eligible to form stars. We compute the probability p» of 
a SPH particle spawning a new star particle of mass m* during a 
time-step At as 



P* = 



M c 



At 



(9) 



where the average star formation timescale t* is obtained by using 
equation ((6}: 



U(P) = 



' M> (P) 

. M 



Myr. 



(10) 



We use m* = m or i g /3 as in IOEFJ051 where m or i g denotes the 
original SPH particle mass. When an SPH particle spawns a star 
particle, the mass in cold clouds is reduced by m*. 

We relax the IRA through a probabilistic treatment of SNe II 
and la. First, we assume that stars more massive than 8Mq explode 
as SNe II. The probability pn of a star particle having an event of 
SN II explosion during a time-step At is given by 



Pn = 



j; +At ru(t')dt' 



t(8M ) 



m(t')dt' 



(11) 



where t is the age of the SSP, r u (t) is the SN II rate for the SSP 
of age t and tfSM©) is the lifetime of the star of original mass 
8Mq. When a star particle satisfies the condition to have an event 
of SN II we distribute mass, metals, and feedback energy, which 
are expelled by the total SNe II between t = and i(8Me) to the 
surrounding gas particles using the SPH kernel weighting. 

Since the SN la rate is much lower than the SN II rate, we 
allow star particles to have only one SN la explosion event for each 
SN la time-step Ati a , which is much longer than the time-steps 
for dynamical calculation. By doing this, we substantially reduce 
the computational expense that is needed mainly for the neighbour 
search to distribute energy, mass, and metals. As for SN II, we can 
compute the SN la probability during a time-step At as 

C +At ri a (t')dt' 

<12) 

where to is the age at which the previous SN la time-step finished 
and ri a (t) is the SN la rate for a SSP of age t. We adopt Ati a = 
100 Myr in this paper. The probabilistic method described above 
statistically relaxes the IRA. 

Finally, we update the mass of the hot phase Mb and the spe- 
cific internal energy of the hot phase u^. The new mass of the hot 
phase Mh' is given by solving the thermal energy equation implic- 
itly using equation Q}, 

A net (p h ',Uh,Z) Mi 



M h ' = M h + AMevp 



-At, 



/Oh' 



(13) 



where AMevp is the mass which evaporates from cold clouds dur- 
ing the time-step, and the new density of the hot phase pi/ is given 
by the mass conservation, Msph = Mh + M c , and by the assump- 
tion that the volume occupied by a SPH particle remains the same 
during the time-step, i.e. 

Msph w , , M h 

= Vol c H , 

P Ph 

where Vol c is the volume occupied by the cold clouds computed 
using the mass-radius relation (eq.[3]l. 

In addition to the adiabatic heating/cooling and shock heating, 
non-hydrodynamic processes such as SN feedback and cloud evap- 
oration also change the specific energy of the hot phase. The new 
specific energy of the hot phase u^' due to the non-hydrodynamic 
processes is also calculated implicitly as 

AQ + (uc - Uh')AM E vp 



«h = Mh + 



Mh' 



(14) 



where AQ is the thermal energy received by the SPH particle dur- 
ing the time-step. 

Note however that we model the SN feedback not as the 
thermal heating of the ISM but as the kinetic heating described 
in the following subsection as winds. This is motivated by the 
results from high-resolu t ion 2D and 3D simulations of the I SM 
( Wad a & Normar] 1200 ll ; IWadal 1200 ll : IWada & Normanl l2007h in 
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which they showed that the ISM is supported by the turbulent mo- 
tion originated in the self-gravity of the gas and galactic rotation. 
They found that e nergy injection b y SNe hardly changes the struc- 
ture of the ISM. iKravtsovl J20030 also indicated that the energy 
feedback from SNe does not change the structure of the ISM but 
significantly reduces the amount of gas in the galaxy in their cos- 
mological simulation performed with an adaptive mesh refinement 
code. 

Lacking the heating by the AQ term, the ISM loses its pres- 
sure quickly and would be fragmented. We therefore introduce min- 
imal heating by imposing the minimum pressure for the star form- 
ing gas as a function of density, 



-Pmin(p) OC p 7off for p> p t h, 



(15) 



by which we mimic the heating from self-gravi t y and galacti c 
she ar motion claimed by IWada & Normarj d200ll) . IWadal hoOlb . 
and Wada & Normanl ( 120070 . We adopt in this paper the value 
7off — 1.4 for the effective adiabatic i ndex, which is close to the 
value found by I Wada & Normanl (2007) for dense gas and is stable 
against gravitational instability (the critical value is 4/3). We nor- 
malise the minimum pressure as P m in(pth) = (7 — l)li4Pth, where 
114 denotes the specific internal energy of the gas with temperature 
10 4 K. Once pressure from the hot phase drops below P m i n , we set 
the temperature of the hot phase to Thot and recalculate the mass 
in cold clouds so that the ambient pressure becomes P m i n . We set 
Thot to 10 6 K in this paper but our results are hardly affected by the 
value of Thot as far as Thot >> 10 4 K. In quiescent disks, most of 
the star forming gas has P m i n by construction. However, for exam- 
ple during a violent merger, the pressure can be significantly higher 
than Pmin due to shock heating and hence the gas can have very 
short star formation timescale. 



2.7 Winds 

As discussed in SH03, we explicitly assume that cold clouds and 
hot ambient medium remain tightly coupled in high-density re- 
gions. Hence, there is no obvious route for the high entropy gas to 
escape from the star forming regions. We thus extend our feedback 
model to accoun t for galactic winds by SN feedback. Our model is 
based on |SH03| and is modified to take the non-instantaneous SN 
explosions into account. 

First, we specify the velocity of the wind, u w , when it leaves 
the disk. Then, we assign the probability p w with which the gas 
particle is added to the wind as 



JJw = 



|A/sphv w 2 
AQ ' 



(16) 



where Msph is the mass of the gas particle and AQ is the feedback 
energy received during the time-step. By doing this, our model be- 
comes insensitive to the numerical resolution, time-step in the sim- 
ulations, and the number of neighbours we use to distribute the 
feedback energy. When a particle is added to the wind, we modify 
the velocity v of the particle according to 

v' = v + n w n. (17) 

For the unit vector n, we choose random orientation along the di- 
rection v x a grav , where a grav is the gravitational acceleration. 
Therefore wind particles are preferentially ejected along the rota- 
tion axis of a spinning object as the 'axial wind' in lSH03l We also 
decouple a new wind particle from hydrodynamic interactions for 
a brief time such that the winds originate from a region close to the 



surface of the star forming region. The full hydrodynamic interac- 
tions are enabled again once the density of the particle has fallen 
to 0.1 pth, or once a time of At = 20 Myr has elapsed, whichever 
happens earlier. 



2.8 Parameter setting and test simulations 

In this paper we use the following parameters for star formation and 
SN feedback as a fiducial model; the threshold density, p t h = 5 x 
10~ 25 g cm -3 , the dimensionless star formation efficiency, c = 
0.005, the dimensionless conductivity efficiency, t^evp = 0.1, the 
SN feedback efficiency, ?7sn = 1, and the wind velocity u w = 500 
km s _1 . The relatively low value of c» compared to other papers 
comes from the fact that we use the dynamical time of each cloud 
not the dynamical time for SPH particles. We choose the values 
of pth and c* so that our model reproduce the observed relation 
between surface gas density and surface star formation rate density 
jKennicutf l998). The most important parameters in our model are 
the feedback efficiency and the wind velocity. T he wind velocity 
we use is comparable to the strong wind model in lNagamine et al.l 
J2004h . For our choice of the parameters and the IMF, the wind 
production rate by SNe II is M w /M* ~ 3. There will be further 
contribution from SNe la at later times. 

In this subsection, we show the behaviour of our model in ide- 
alised simulations of disk formation from virialised gas in the ex- 
ternal NFW potential ( tNavarro et al J 19960 . Here, M vlr is the virial 
mass of the system and we always assume that 10 per cent of the 
mass is in baryonic form. The initial angular momentum in terms 
of the spin parameter, A = J|P| 1/2 /(GM 2 / r 5 ), j s se t to 0.07 with 
the assumption that the specific angular momenta j(r) of spherical 
shells are all aligned and their magnitude is given by 



j(r) oc 



M(r) 

Mvir 



In Fig. [T] we show the star formation rates in the simulations 
for M v i r = 10 12 Mq with three different resolutions: iV gas = 
10000, 50000, and 250000. The star formation rates quickly con- 
verge. The simulation with -/V gas = 50000, which is a compara- 
ble resolution that we will use in cosmological simulations, has 
very similar star formation history to that in the simulation with 
N gas = 250000. Hereafter we thus only show the simulations with 
iV gas = 50000. 

Now we show how star formation rates depend on the values 
of feedback parameters, 7/sn and v w . We here run the simulations 
with two additional models: the slow wind model (h]sn = 1 and 
w w = 250 km s~ ) and the weak feedback model (rysN = 0.5 and 
w w = 250 km s" 1 ). The wind speed in the slow wind model is half 
of that in the fiducial model, therefore it has an uncomfortably high 
mass loading factor r] w ~ 12. Of course the mass loading in the 
weak feedback model becomes half. 

Fig. H] shows the star formation rates in galaxies with these 
three wind models for halos of M vil = 1O 12 M and 1O 1O M . 
Since w w = 500 km s _1 exceeds the escape velocity of the halo 
with M vir = 10 12 M Q and u w = 250 km s -1 does not, the star for- 
mation is more strongly suppressed in the fiducial model in spite of 
the same feedback efficiency. We find that the star formation rates 
of the slow wind and weak feedback models have regularly time 
spaced peaks of activity. For the adopted NFW potential, the wind 
particles with the initial velocity, v w — 250 km s _1 , will return to 
the disk in ~ 0.13 Gyr. This time-scale might be responsible for 
the periodicity seen in the star formation rates. The star formation 
rate is highest in the weak feedback model for this halo as expected. 
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Figure 1. Star formation rates in simulations of disks growing within ha- 
los of M v i r = 1O 12 M0, with different resolutions. The simulations with 



10000, 50000, and 250000 are indicated by the solid, dotted, and 



dashed lines, respectively. 



However situation changes for the halo of M v \ r — 10 10 Mq. Now 
the mass loading factor is more important than the wind velocity 
because both wind velocities exceeds the escape velocity of this 
halo. Hence the star formation rate is highest in the fiducial model 
and lowest in the slow wind model although all models have star 
formation rates ~ 10~ 3 Mq yx -1 after 1 Gyr. Since we have to 
suppress star form ation even in relatively large halos to form disk 
galaxies (OEFJ05), we employ the fiducial parameter set in our cos- 
mological simulations. 

In Fig. [3] we show the relations between the star formation 
rate per unit area and the surface gas density for the fiducial (top 
panel) and weak feedback (middle panel) models. We estimate the 
surface star formation rate densities by computing the surface den- 
sity of stars younger than 3 x 10 7 yr in cylindrical bins. Both 
models show rea sonably good agreement with the target relatio n 
jKennicuttll 19981 ; the target relation we show is eq. (25) of lSH03l) . 
The use of different wind parameters do not affect this relation in 
our simulations, while the ISM can have higher gas density in the 
weak feedback model. For the galaxy in M v i r = 10 10 Mq halo, 
winds are so effective that there are only a handful of young stars 
in the galaxy. We hence only show, in the bottom panel of Fig. [3] 
the relation between the global surface gas density and the global 
surface star formation rate , which are averaged over the galaxy. 
The different points in this panel correspond to the relations at dif- 
ferent epochs. We find that the global relation in the 10 10 Mq halo 
is also in reasonable agreement with the target relation. 



3 BLACK HOLE ACCRETION AND AGN FEEDBACK 
3.1 Seed black holes 

In the cosmological sim ulations we run a friends-of-friends halo- 
finder jPavis etal.ll 19851) 200 times on the fly from z = 15 to 0. 
When we find a halo that contains more than iV t h dark matter parti- 
cles and no BH particles, we put a seed BH on a density peak of the 
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Figure 2. Star formation rates in simulations of disks growing within halos 
ofM vir = 10 12 M Q (upper panel) and M vir = 10 10 M Q (lower panel). 
The solid, dotted, and dashed lines indicate the fiducial, slow wind, and 
weak feedback models, respectively. 



stellar distribution in the halo. We assume 260 Mq as the mass of 
the seed black hole. This mass is consistent with the theory of stel- 
lar evolution, which shows that the nuclear burning in very massive 
stars above 260 M q is unable to halt gravitational collapse (e.g. 
iHeger et alj|2003» . Note that as far as the seed BH mass is suffi- 
ciently small, our results hardly change because we do not impose 
an upper limit on the mass accretion rate. Throughout this paper 
we use Nth = 128. The simulation employing Nth = 256 yields 
almost the same results except for the total number of BHs. The 
insensitivity on the choice of these parameters relies on our AGN 
feedback model that is only effective when the mass accretion rate 
onto the BH is much lower than the Eddington rate (see i]3.3t . 



3.2 Mass accretion by the radiation drag 

One mechanism yielding the proportionality of the BH mass 

to its host bulge mass has been propos ed by lUmemural 

fcOOlh and Kawakatu & Umemural d2002l) . According to 
iKawakatu & Umemural in the central regions of galaxies the drag 
due to stellar radiation may result in a loss of angular momentum 
of a clumpy ISM and gives rise to mass inflow toward the centre. 
The optical depth of a gas cloud is f = XdPcTc — Xd m c/r 2 , 
where p c , m a , and r c are the density, mass, and a size of a cloud. 
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Figure 3. Star formation rate per unit area versus gas surface density in an 
idealised simulation of disk formation in a halo with M V \ T = 1O 12 M0 (top 
and middle panels). The surface star formation rate densities are computed 
using the surface density of stars younger than 3 X 10 7 yr in cylindrical 
bins. The symbols represent star formation rates at t = 0.5, 1, 2, 3 Gyr, 
indicated by plus signs, diamonds, triangles, and squares respectively. The 
top and middle panels respectively show the fiducial and weak feedback 
models. The galaxy in the halo with Af v ; r = 10 10 Mq is shown in the 
bottom panel, and there both the surface star formation rate density and 
the surface gas density are averaged over the galaxy. Points represent the 
evolution of the relation from 0.25 Gyr to 3 Gyr with a time-step of 0.25 
Gyr. The solid and dotted lines indicate the target relation and the cut-off 
surface density respectively. 



In principle we can use the cloud model we described in §272] We 
however assume that the clouds are all identical and randomly 
distributed in a star forming region as in iKawakatu & Umemural 
for simplicity. The mass extinction coefficient Xd is given by 
Xd = 300 cm 2 g- 1 (a d /0.1/an- :L )(p B /g cirr 3 )(Z/O.3Z ), 
where a<j is the grain radius, p s is the density of solid material 
within the grain (e.g. ISpitzeijll97ol) . and Z is the metallicity of 
gas. We use a<j = 0.1 /im, p B = 1 g cm -3 , and Z is directly 
obtained in the simulations as the mean metallicity of the gas in 
a star forming region. The total optical depth of a star forming 
region can then be approximated as 



tsfr(£) 



3Xd M c (t) 

4tt RsFR(t) 2 



(18) 



where Rsfr and M c are the size of a star forming region and the 
total mass of the clouds in the star forming region. 

The mass accretion rate due to the radiation drag mechanism 
described above is well approximated by 



where Lsfr is the total luminosity of stars in the star form- 
ing region. We calculate Lsfr by summing up the bolometric 
luminosities of all stars, each of which is obtained as a func- 
tion of its age from a look-up table generated by PEGASE2 
jFioc & Rocca-Volmerangei 1 1997b . We here ignore the metallic- 
ity dependence of the stellar luminosity and use a table for the 
solar metallicity becaus e the dependence on metallicity is weak. 
IKawakatu & Umemural found that the efficiency r^drag is maxi- 
mally 0.34 in the optically thick regime assuming a spherical star 
forming region with a uniform density. In principle ?7drag can be 
much larger or smaller depending on the geometry and density 
structure of the star forming region (Kawakatu private communica- 
tion). We thus change the efficiency and adopt 7/drag = 1 in this pa- 
per, considering more centrally concentra ted density distribution i n 
real star forming regions than assumed in Kawakatu & Umemural 

Once we define the radius of a star forming region around a 
BH, iisFR, we obtain \d, M c , and Lsfr in equations d 1 8fc and dl9l l 
directly in simulations, and hence we can compute tsfr and Mdrag 
self-consistently; the effects of outflows are also implicitly incorpo- 
rated in the estimate of these variables because we calculate them 
self-consistently at each time-step of the simulation. In order to de- 
fine the size of the star forming region, -Rsfr we start from a sphere 
around a BH that contains 40 SPH particles, and then we search the 
radius of the sphere that maximises Mdrag by increasing the radius 
of the sphere. We employ this radius as the size of the star forming 
region. When the mean gas density of the initial sphere is lower 
than the threshold density for star formation, p t h, we consider that 
there is not enough ga s around the BH to fue l accre tion and we 
simply set A/drag = 0. Kawak atu & Umemural {2004) argued that 
radiation from bulge stars contributes to the mass accretion to the 
central hole more efficiently than that from disc stars. We crudely 
mimic this effect by using spherically averaged values to define star 
forming region; by doing this we put less weight to the star forming 
gas that has a flattened distribution. We will show how star forma- 
tion in the disc and bulge correlate with the mass accretion rate, 
Mdrag. in our simulations later. 

Strictly speaking, one should distinguish Md rag from the ac- 
cretion rate onto the BH, M, because the radiation drag is not likely 
to remove the angular momentum thoroughly; some residual an- 
gular mom entum will terminate the radial contra ction of the ac- 
creted gas dSato et alj|2004l) . iGranato et al.l d2004l) introduced the 
viscous timescale in which the mass accumulated at the galactic 
centre by the radiation drag is fed to the BH. We however assume 
M — ALjrag for simplicity. The time delay in mass fuelling intro- 
duced by the viscous timescale may become important when we 
compare the luminosity of AGNs and their host galaxies in detail. 



3.3 AGN feedback 



Motivated 



( Croton et 



the success of rec ent semi-analytic models 
200o1 ; iBower et ai1l2006l) . we assume that AGNs di- 



rectly heat diffuse hot gas confined in dark halos through the pro- 
duction of jets. As the jet produ ction mechanism, we employ the 
model proposed by Meier (2001 ) which is an hybrid version of the 
Blandf ord-Znaiek jBlandford & Znaiekll 19771) and the Blandford- 
Payne (Blandfo rd & Pavnd Il982h processes. In this hybrid model 
the jet power is generated by the magnetic field threading the ac- 
cretion flow inside and outside the the ergosphere of the BH. This 
model is able to draw upon the rotational energy of the accretion 
flow as well as the spinning black hole in order to drive jets, and its 
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basic features are supporte d by recent relativistic nu merical simu- 



lations of jet form ation (e.g. Hawlev & Krolik 2006b). We refer the 
read ers to[M eier ( 200 1 ) and Nem men et al J J2007t) for more details. 

iMeierl ( 1200 II) presented models both for non-spinning 
(Schwarzschild) and for spinning (Kerr) BHs. Semi-analytic cos- 
mological simulations of the spin evolution of BHs through merg- 
ers and gas accretion dVolonteri et al.l 2005) and estimates of the 
radiative efficiencies of glo bal populations of quasars based on 
Soltan-type argu ments (e.g. ISoltan | 19821 ; | Yu & Tremainel [2002 ; 
IWang et all 120061) sugge st that most nearby massive BHs are 
rapidly rotating. Recently INemmen et al.l J2007) also reported that 
BHs in nearby elliptical galaxies are likely to have high spins 
j ~ 0.7 — 1. W e thus only consider spinning BHs in this pa- 
per. lMeierld200ll) parameterised the jet luminosity for standard thin 
discs as 



1n 42.7 

10 erg s 
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and for RIAFs as 



r riaf _ in45.i -i/o?riaf\ / rn 

jet w 10 ergs l^^J m9 lai 



x (0.55/ 2 + 1.5/j+j 2 ), 



(20) 



(21) 



where qsd and qriaf are the viscosity parameters for standard 
thin discs and for RIAFs, respectively, mg is the black hole mass 
in units of 10 9 Mq, rh = M / Msaa is the accretion rate scaled 
to the Eddington rate assuming the 10% energy conversion effi- 
ciency (A/sdci = 22 trig Mq yr _1 ), and / and g are the ratios 
of the actual angul ar velocity and t he azi muthal magnetic field to 
th ose calculated by iNaravan & Yil j 19951) respectively. According 
to lMeierl d200ll) we employ / = 1 and g = 2.3. Since the limited 
numerical resolution of our cosmological simulations prevents us 
from resolving the details of the propagation of jets and its impact 
on the surrounding medium, we simply assume that a fraction of 
the jet energy A_Bp B = r/AGN-t/jct Ai is delivered as thermal en- 
ergy to the neighbouring 40 diffuse gas particles that have density, 
p < O.lpth- In our simulations, the size of the region heated by 
the AGN feedback is typically ~ 10 kpc. Unless otherwise stated, 
we use cisd = oriaf = 0.1, j = 0.5, and ^agn = 0.1. The 
adopted value of j corresponds to the intermediate case between 
Schwarzschild (no-spin) and maximally rotating BHs. 

As previously mentioned, we adopt two distinct regimes of 
accretion flows: standard thin discs (optically thick, geometrically 
thin, radiatively efficient) and RIAFs (optically thin, geometrically 
thick, radiatively inefficient). The parameter which controls the 
state of the accretion flow is rh, and the critical value of the accre- 
tion rate which sets th e division between t he two regimes is given 
by m cr it ~ a 2 (e.g. INaravan et alJI 19980 . Since the RIAF solu- 
tion ceases to be valid above rhcrit, for rh ^ rh cr it a RIAF exists 
and for rh > rh cr it a thin disc occurs. As a consequence of the 
adopted theoretical model of jet production, we have strong AGN 
feedback only for RIAFs and the maximum AGN feedback occurs 
at rh = rhcrit ■ For the parameters we chose, we have the following 
jet production efficiency in RIAFs, 



LjJ|; AF ~ 2.6 x 10 1 M<? for rh rh cr i t , 

which is much higher than that in standard discs, given by 



8.1 x 10 Mc for rh > rhcrit, 



(22) 



(23) 



and J21b correspond to the "radio-loud mode" in our simulations, 
and the equations J23b and i20l represent the "radio-quiet mode". 

Note that, in principle, we could have even higher efficiencies 
of jet production, because the energy is drawn from the "spin en- 
ergy" stored in the BH. We explore the case of a higher efficiency 
corresponding to a higher BH spin later in <j5 .3 1 



4 COSMOLOGICAL SIMULATIONS OF DISC GALAXY 
FORMATION WITH AGN FEEDBACK 

To study the role of AGN feedback in disk gala xy formation, we 
use the initial conditions presented in IOEFJ051 from which they 
produced an extended disk galaxy. The background cosmology is 
the so-called ACDM with the mean matter density Ho = 0.3, Hub- 
ble parameter, h = .ffo/100 km s _1 Mpc _1 = 0.7, cosmological 
constant term, ^Ia = Ao/(3Hq) = 0.7, amplitude of mass fluc- 
tuations, as = 0.9, and mean baryon density, fib = 0.04. The 
size of the periodic simulation box is 35.325h~ 1 Mpc and the ini- 
tial redshift is z — 49. We put the high-resolution dark matter 
particles and gas particles in the region where a halo with mass 
M vir = 1.2 x lO 12 h _1 Af forms at z = 0. The region external 
to this is populated with high-mass dark matter particles (bound- 
ary particles), the function of which is to reproduce the appro- 
priate tidal fields. The circular velocity, spin parameter and col- 
lapse redshift of the selected halo are « c (r v ir) = 155 km s _1 , 
A = J\E\ 1/2 /(GM 5/2 ) = 0.038 and z c ~ 1.5, respectively, 
where z c is defined as the redshift at which the main progenitor 
had half the final halo mass. 

The masses of the SPH and high-resolution dark matter parti- 
cles are ~ 2.6 x 10 6 and ~ 1.7 x 10 7 h~ 1 M Q , respectively. Thus 
the N th = 128 implies that halos larger than ~ 2.6 x lO 9 /i _1 M 
are allowed to have BHs. The gravitational softening length are kept 
fixed in comoving coordinates for z > 3; thereafter they are frozen 
in physical units at a value (of the equivalent Plummer softening) of 
e = 0.5 and 1 kpc for the baryonic particles (SPH particles, stars, 
and BHs) and high-resolution dark matter particles, respectively. 
The gravitational force obeys the exact r -2 law at r > 2.8e. 



4.1 Black hole particles 

When we find a halo with Ndm ^ Nth that does not contain any 
BH particles and is not contaminated by the boundary particles, we 
turn the star particle that has the largest stellar density into a seed 
BH. We thus discriminate the BH mass from the mass of the BH 
particle (particle mass). The former is updated every time-step ac- 
cording to the mass accretion rate obtained from the radiation drag 
model. Once BH mass exceeds its particle mass, the BH particle 
increases its particle mass by absorbing gas particles n in the star- 
forming region around it with the probability: 

M BH - M p 



1 



iVsFR M n ' 



(24) 



where in equation {23} we adopt mg = rh — 1. The equations d22t 



where AIbh and Af p are, respectively, the BH mass and the mass 
of the BH particle, A^sfr is the number of gas particles in the star- 
forming region, and M n is the mass of the gas particle n. 

BH particles can increase their BH masses and particle masses 
by BH mergers. However, how long it would take to harden a BH 
binary until finall y gravitational wave emission becomes i mpor- 
tant is still unclear dMakino & Funatol2o64l : lEscala et alj2004D . We 
simply assume that BHs can merge into a single BH if two or more 
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BHs are within a gravitational softening radius and they are grav- 
itationally bound. Therefore the BH merger rates obtained in our 
simulations are likely to be the maximum possible rate. Techni- 
cally, we apply the friends-of-friends algorithm to find groups of 
BHs which may potentially merge originating new BHs. 

An interesting problem arises owing to the finite numerical 
resolution. In reality, supermassive BHs are much heavier than stars 
and dark matter. Therefore they are quickly brought to the galac- 
tic centres and stay at the bottom of the potential wells owing to 
the dynamical friction. In cosmological simulations, the mass of a 
BH particle is comparable to those of gas and stellar particles and 
lighter than the dark matter particles in the early stage of its evolu- 
tion. As a result, BH particles tend to oscillate around galactic cen- 
tres with the velocity dispersions of the stellar systems. To circum- 
vent this problem, we allow BH particles to behave as peak tracers 
at galactic centres. At each time-step, we compute local stellar den- 
sity fields within the gravitational softening radii of BH particles, 
e. We then displace the BH particles by a small distance, Al, along 
the local density gradients. After several tests, we chose 

Al = MIN(0.01e, 0.03 \v\ At), (25) 

where v is the velocity of a BH particle and At is the time-step. 
This Al correction sufficiently suppresses the oscillation of the BH 
particles at galactic centres. 

We adopt a maximum time step for BHs, Ai max = I Myr, 
in order to follow mass accretion and associated feedback with a 
reasonable time resolution. Hence a time step for a BH particle is 
defined as A£bh = min(Af max , Aid yn ), where Aidyn is a time 
step determined by the dynamics. 

4.2 Models 

We first compare three simulations in order to show the effects of 
AGN feedback. The first simulation does not have feedback from 
AGN. We refer this model as the 'no-fb' model. In the second 
model, we distribute 1% of jet energy to the ambient halo gas par- 
ticles defined by their densities, p < 0.1p t h in the form of thermal 
energy using SPH smoothing. Since only 1% of jet energy is avail- 
able for heating, we call this model the 'weak-fb' model. In the 
'reference' model, we adopt the same AGN feedback scheme as 
in the weak-fb model but assume 10% of jet energy is available to 
heat the ambient halo gas. We will explore other models in which 
we change the implementation of AGN feedback and strength of 
SN feedback in $53\ 



5 RESULTS 

To give a visual impression of our simulations, we show some 
snapshots of the simulation with AGN-radio feedback ('reference' 
model) in Fig. [4] There are a number of merger events at high red- 
shift (z > 1). The disk formation starts at z ~ 1, and finally a disk 
galaxy forms. A lot of BHs form in the small halos and they are 
brought to the main progenitor by merging satellites. The BH parti- 
cles nicely trace the centres of the stellar objects resulting from our 
correction applied to the positions of BH particles. There are also 
many BHs which are not associated to any stellar objects. They 
are scattered into the halo when small satellites, which harbour a 
small BH at their centres, are destroyed during mergers. Most of 
the drifting BHs hence have very small mass. 

In the following sections, we present the main results of the 
paper. In §5.1, we compare the results of simulations run with 




Figure 5. Projected stellar distributions within 50 h~ 1 kpc boxes centred 
on the galaxies at z = 0. The left, middle, and right panels show the 'no-fb' , 
'weak-fb', and 'reference' simulations, respectively. The viewing angles are 
chosen to be edge-on for upper panels and face-on for lower panels. The 
brightness indicates the projected stellar mass density, and the same scaling 
is used for each simulation. 

the standard feedback parameters (the 'reference simulation') with 
simulations using weaker AGN feedback (the 'weak-fb' model) 
and no AGN feedback at all (the 'no-fb' model). In §5.2 we look at 
the importance of the AGN feedback for the bolometric luminosity 
of the forming galaxy, while in §5.3, we explore the dependence 
of the results on the way in which the feedback is coupled to the 
surrounding gas and the impact of assuming a lower efficiency of 
the feedback from SNe. 

5.1 Simulations with different AGN feedback efficiencies 

Now we compare galaxies formed in three simulations, i.e. the 'no- 
fb', 'weak-fb', and 'reference' simulations. We first show the stellar 
distributions within 50/i _1 kpc boxes centred on the galaxies at 
z = in Fig. [5] The edge-on and face-on projections are selected 
to be perpendicular and parallel to the angular momentum vector of 
the stellar component within the central 10ft _1 kpc sphere. In all 
of the models, the galaxy has an extended stellar disk. The 'no-fb' 
and 'weak-fb' models are very similar. From the face-on view, it 
is found that the 'reference' galaxy has a disk with lower surface 
stellar density. 

Fig. |6] shows the I band surface brightness profiles for these 
galaxies. We also show the surface brightness profile of a disk 
galaxy obtained by OEFI05 using the same initial condition. An in- 
ter esting dif ference from the surface brightness profile of the galaxy 
by|QEFJ05 can be seen at the centre. Our galaxies do not show the 
extreme central concentration. The origin of this difference is the 
different implementation of SN feedback. While we implement it 
as kinetic winds, they deposited feedback energy in the form of 
thermal energy. Therefore they had to calculate the pressure gradi- 
ents caused by the thermal feedback to blow the gas away. Unfortu- 
nately, the SPH does not allow steep gradients in physical quantities 
by definition unless one uses a sufficiently high resolution. 

All our three galaxies show similar central surface brightness. 
The main difference due to AGN feedback can be seen in the sur- 
face brightness of the disks. While 'no-fb' and 'weak-fb' galaxies 
show the almost same surface brightness profiles, the 'reference' 
model has a less extended disk. We do not try to determine the 
bulge-disk decomposition based on the surface brightness profiles 
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Figure 4. Snapshots of the 'reference' model (x—y projections). The rows correspond to 2 = 6, 3, 2, 1, and from top to bottom. The first and the second 
columns show gas and stellar distributions in the comoving 3 ft -1 Mpc cubes centred on the main progenitor, respectively, where the brightness represents 
three dimensional density. The third and the fourth columns show gas and stellar distributions in the physical 50 hT 1 kpc cubes, respectively. The circles in 
each panel indicate the position of the BH particles and their radii are proportional to the log of the BH masses. 
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Figure 6. Surface brightness profiles in the / band for the 'no-fb' (pluses), 
'weak-fb' (asterisks), 'reference' (diamonds) galaxies, and for the galaxy 
with the 'shock-burst' model of Okamoto et al. (2005; dotted-line). The 
solid lines are exponential fits for r < 15 kpc with the scale lengths 5.4, 
5.6, and 4.0 kpc for the 'no-fb', 'weak-fb', and 'reference' galaxies, respec- 
tively. 
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Figure 7. The formation history (star formation rate as a function of time) 
of the stars that lie within 25h~ 1 kpc from the galactic centre at z = 0. The 
solid and dotted lines correspond to the galaxies in the simulation without 
and with BHs, respectively. 



since they provide only partial information about the morphology 
of galaxies. Belo w, we will show th e decomposition based on the 
dynamics of stars dAbadi et alj|2003r) . This is a more accurate way 
to characterise the relative importance of the bulge and disk com- 
ponents. 

To investigate effects of AGN feedback quantitatively, we 
present the star formation histories of the stars that lie within 25/i _1 
kpc from the galactic centres at z — in Fig. [7] All galaxies have 
qualitatively similar star formation histories and have peak star for- 
mation rates at z ~ 1.5. During the starburst {t = 2.5 — 6 Gyr), the 
'reference' galaxy shows lower peak star formation rate compared 
with other galaxies. As we will show later, the AGN is not a RIAF 
during the starburst. Hence AGN feedback from a RIAF is not re- 
sponsible for this lower star formation rate. The feedback from 
standard disk (equation d23b ) might suppress the starburst. However 
we should not overinterpret the differences during the starburst be- 
cause the system is strongly non-linear during this epoch and hence 
small changes can result in relatively large difference. In fact, the 
model which employs stronger AGN feedback has a higher star for- 
mation rate during the starburst than in the 'reference' model as we 
will show later in Fig.l 151 

After the starburst, the 'reference' galaxy has significantly 
lower star formation rate than other galaxies. This suppression of 
star formation is likely to be caused by AGN feedback from a RIAF. 
The lower star formation rate in the 'reference' galaxy explains its 
less extended disk in Fig. [6] because it is the star formation after 
z ~ 1 that builds up the stellar disk (see Fig. |4). The 'weak-fb' 
galaxy has almost the same star formation rate as the 'no-fb' galaxy. 
Hence it can be said that, in order to give a visible impact on galaxy 
formation, at least ~ 10% of jet energy has to be given to the halo 



gas. We will later explore the case in which the jet energy is not 
given to the diffuse halo gas but delivered to the surrounding ISM. 

Since the strength of AGN feedback in our model is closely re- 
lated to the central BH mass and the accretion rate onto the BH, we 
plot the evolution of the central BH of the main progenitor and the 
mass accretion rate in units of the Eddington rate in Fig.[8]in each 
model. Note that the mass accretion rate M is different from the 
mass growth rate, A/bh ; the latter includes the mass increase due to 
BH mergers. In all simulations, the central BHs rapidly evolve un- 
til t ~ 5.5 Gyr at which the starburst ceases. The final BH masses 
reach ~ 10 7 Mq. The stepwise jumps seen in the evolution of the 
BH mass indicate the mass growth by mergers. 

The consumption of the available gas for star formation ex- 
plains the decline of the accretion rate onto the central BH. At 
t ~ 5.5 Gyr m reaches the critical accretion rate m cr it, below 
which the accretion disk becomes a RIAF. Our model has the max- 
imum AGN feedback when the mass accretion rate is equal to the 
critical accretion rate. This epoch is indicated by the crossing of 
the solid line and the horizontal dashed line in the lower panel of 
Fig. [8] Once the accretion disk becomes a RIAF after z < 1 , in 
the 'reference' model, the accretion rate drops to zero rapidly and 
there is no mass accretion after the strongest production of AGN 
feedback at t ~ 5.5 Gyr except for some episodic accretion events. 
This implies that there is almost no star formation activity at the 
centre of the main progenitor at z < 1 for this galaxy. On the other 
hand, lacking the strong AGN feedback, there are not such sudden 
decrease in accretion rates in the 'no-fb' and 'weak-fb' models. 
The gradual decline is caused by consumption of gas due to star 
formation and associated SN feedback( winds). 

Our implementation of AGN feedback cannot directly sup- 
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Figure 8. Upper panel: the mass growth of the central BH of the main progenitor. Lower panel: the mass accretion rate onto the central BH of the main 
progenitor in units of the Eddington rate (m = M / M^dd)- The horizontal dashed line indicate the critical mass accretion rate m cr ; t below which the 
accretion flow becomes radiatively inefficient (RIAF). The 'no-fb', 'weak-fb', and 'reference' simulations are shown from left to right panels. 



press the star formation because it injects energy not into the star 
forming dense gas but into the diffuse halo gas (p < 0.1/9 t h)- Thus 
it is the suppression of the cooling of the halo gas around the centre 
of the galaxy which is responsible for the halting of star formation 
around the BH in the 'reference' simulation. The impulsive feature 
seen in the accretion rate implies the existence of a self-regulated 
cycle that operates between gas cooling, star formation, and pro- 
duction of jets in the RIAF. The fact that the central surface bright- 
ness is almost the same in the three models in spite that the galaxy 
with AGN feedback hardly form stars at its centre after z ~ 1 sug- 
gests that the vast majority of stars at the centre have already been 
formed at z ~ 1 in all galaxies. 

We now return to the separation of the disk and bulge compo- 
nents. As noted earlier, it is best to separate the bulge and disk com- 
ponents using a dynamical decomposition ( Abadi et ah 20031). For 
this, we first compute the angular momentum J z of each star par- 
ticle parallel to the net angular momentum of stars within 10/i -1 
kpc, and the angular momentum of the corotating circular orbit, 



J C (E), where E is the energy of each particle. The ratio J z /J c 
defines an orbital circularity. In Fig. [9] we show the probability 
distribution of this orbital circularity for stars within 25/i _1 kpc 
from the centre of each galaxy. A disk component should have 
Jz/Jc(E) ~ 1 and we see that such a component clearly d om- 
inates in each galaxy. Comparing to the results by OEFJ051 the 
suppression of the formation of spheroidal comp onents seems to 
rely la rgely on the strong SN feedback by winds. rGovernato et al.l 
(2007) also produced disk dominated galaxies in halos with quieter 
merger histories (in which the epoch of the last major merger is 
zimm > 2) using a different feedback recipe from ours. To iden- 
tify the stars in the spheroid, we assume a non-rotating spheroid, 
i.e. that stars in the spheroid are symmetrically distributed around 
zero in Fig. [9] where all stars having J z / J C (E) ^ are identified 
as a half of the spheroidal component. Another half of the spheroid 
with Jz/J C (E) > are defined statistically so that the total angular 
momentum of the spheroid becomes zero. All remaining stars are 
identified as the disk component. As expected from the similarity in 
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Figure 9. Mass-weighted probability distributions of the orbital circular- 
ity, J z I J C (E), of stars within 25/i — 1 kpc from the galactic centres. The 
solid, dotted, and dashed lines indicate the 'no-fb', 'weak-fb', and 'refer- 
ence' galaxies, respectively. 



the star formation histories, all galaxies have similar distributions 
in this plane while the disc component in the 'reference' model is 
slightly less significant because of the relatively strong AGN feed- 
back. 

In Table Q] we show the total and spheroid's masses, total lu- 
minosity, and the disk-to-total ratios in mass, g band, and r band. 
Since colour of a stellar population is sensitive to its metallicity, we 
here include the metallicity dependenc e to compute the luminosity 
at eac h passband by using PEGASE2 dFioc & Rocca-Volmerangel 
1 1997b . The masses of the spheroidal components directly reflect 
amplitude of starbursts and therefore the 'no-fb' galaxy has the 
most massive spheroidal component. Nevertheless, all galaxies 
have disk-dominated morphology. 

It is interesting to see when each component forms and how 
the formation history of each component is related to the accretion 
history of the central BH. In Fig.[l0] we show star formation histo- 
ries of the disk and spheroidal components in each galaxy at z = 0. 
We also show the rescaled accretion rate onto the central BH of the 
main progenitor of each galaxy. Note that the star formation rates 
include star formation in all building blocks (all progenitors) but 
the mass accretion rates are only the rates onto the central BHs of 
the main progenitors. Thus at high redshift, the mass accretion rates 
presented in the figure considerably underestimate the net accretion 
rates. Since the galaxies do not undergo any significant mergers af- 
ter z ~ 1, it is safe to compare the star formation rates and the mass 
accretion rates at low redshift. We find that there are good correla- 
tions between the star formation rates of the spheroid and the mass 
accretion rates onto the central BHs at low redshift, in particular 
between t ~ 5.5 and 8 Gyr. After the starburst, the star forma- 
tion rates of the spheroids rapidly decreases in all simulations. The 
accretion rates onto the BHs also decrease at the same time. The 
contribution from disk stars to the mass accretion rate due to the 
radiation drag effect should be largest during the peak of star for- 
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Figure 10. Star formation histories for the disk and spheroid stars within 
25h~ 1 kpc from the centre of each galaxy at z = 0. The 'no-fb', 'weak-fb', 
and 'reference' galaxies are shown in the top, middle, and bottom panels, 
respectively. The star formation histories for the disk components are repre- 
sented by the solid lines and those for the spheroidal components are indi- 
cated by the dotted lines. We also show the mass accretion rates, which are 
multiplied by 1000, onto the central BHs in the main progenitors (dashed 
lines). Note that while star formation not only in the main progenitors but 
also in all other building blocks is counted in the star formation histories, 
the mass accretion rates only represent the rates onto the central BHs of the 
main progenitors. 

mation in the disc component (t ~ 5 Gyr) because at this epoch 
the disc is gas rich and very compact. We have however confirmed 
that even at this epoch the contribution from disc stars is less than 
20% and most of the time less than 5% in all models. Considering 
the fact that half of the spheroid stars are defined statistically, the 
correlation between the star formation rate and the mass accretion 
rate of the spheroid in each galaxy is significant. After t ~ 8 Gyr 
there are many short lived star formation episodes in the spheroids, 
which do not correspond to gas accretion events. These star forma- 
tion events occur in small satellites before they become a part of 
the spheroidal components by minor mergers. 

Next we show in Fig. [TT] the evolution of the relations be- 
tween mass of the central BHs of the main progenitors and the 
mass of the host spheroids at z = 3, 2, 1.5, 1, 0.5, 0.2, and 0. 
We here define the main progenitors as stars within 10% of the 
virial radii at each redshift. We have confirmed that the spheres 
with radii of 0. lr v i r contain the stellar systems defined as main pro- 
genitors by the friends-of-friends algorithm in S]5.2I We have also 
checked by eye that at selected redshifts the main progenitors are 
not undergoing major mergers and the spheres do not contain mas- 
sive satellites, both of which ruin the dynamical decomposition. 
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Table 1. Total and spheroid's mass, g and r band total luminosities, and disk-to-total mass and luminosity ratios for the galaxies at z 



= 0. 



A/tot/ M© Msp/M @ M g M r (§) 



no-fb 6.01 X 10 10 1.98 X 10 10 -20.8 -21.3 0.67 0.83 0.79 

weak-fb 5.98 x 10 10 1.44 x 10 10 -20.9 -21.3 0.76 0.88 0.85 

reference 4.32 X 10 10 1.21 X 10 10 -20.6 -21.0 0.72 0.86 0.83 




Figure 11. The relation between mass of the central BHs of the main pro- 
genitors and the stellar mass of the host spheroids at z = 3, 2, 1.5, 1, 0.5, 
0.2, and 0. The relations for the 'no-fb', 'weak-fb', and 'reference' simula- 
tion are shown by the plus signs, diamonds, and triangles respectively. The 
solid line indicates Aip,H = O.OOlAfgp. 



We then compute the mass of the spheroidal component of each 
main progenitor based on the dynamical decomposition, noting that 
this method of decomposition is less vulnerable to insufficient res- 
olution compared with the decomposition based on surface density 
profiles. This feature is important at high redshift where the size of 
the galaxies are much smaller than that at z = 0. Fig. QT] shows 
that the Mbh-Msp relations evolve similarly in all simulations be- 
cause the BHs and spheroids mainly increase their mass during the 
starbursts and our AGN feedback is not effective during this phase. 
Before the starburst, the BH mass is far below the local relation, 
Mbh ~ 10 _3 Msp, and then it rapidly moves towards the local 
relation during the starburst (z — 2 ~ 1). The delay in the mass 
evolution of the BHs compared to the spheroids is resulting from 
the fact that the mass accretion rate due to the radiation drag model 
depends on the metallicity of the ISM. The mass accretion by this 
model is not efficient until the ISM is metal enriched and the opti- 
cal depth becomes ~ 1. After the starburst, the mass of BHs stays 
almost constant as seen in Fig.[8] The relation is imprinted in the ra- 
diation drag model because the model provides a good correlation 
between star formation in a spheroid and the mass accretion onto 
a BH. It is important to note that the radiation drag model predicts 
considerably lower mass ratio between central BHs and their host 
spheroids before the starburst. 
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Figure 12. Colour-magnitude diagrams for the 'no-fb' (left), 'weak- 
fb' (middle), and 'reference' (right) galaxies. The galaxies at z = 
3, 2, 1.5, 1, 0.5, 0.2 and are indicated by the plus signs connected by 
the dotted line. The colours of their spheroidal components are also indi- 
cated by the diamonds. The upper and lower thin solid lines represent mean 
colours of the red and blue populations taken from Baldry et al. (2004). 
The u — r colour of simulated galaxies are converted to the SDSS colours 
by using (it - r)(SDSS) = (u - r)(AB) + 0.05 (Abazajian et al. 2003) to 
compare with the data by Baldry et al. (2004). 



In order to illustrate how the colours of our galaxies evolve, 
we show the evolutions of the main progenitor of each galaxy 
and its spheroid in a colour-magnitude diagram (Fig. 1 1 2b - The 
mea n colours of the re d and blue populations of the SDSS galax- 
ies JBaldrv et al J 1 20041) are indicated by solid lines. All galaxies 
show similar evolution on the colour-magnitude plane. This con- 
firms that the AGN feedback has only minor effects on evolution 
of disk galaxies. The colour of the galaxies and their spheroids at 
2 = are consistent with the observed colours of local blue and 
red galaxies, respectively. 



5.2 Relative importance of the AGN 

In order to show the relative importance of the AGN in galaxy for- 
mation, we show the bolometric luminosities generated by stars in 
the main progenitors and the bolometric luminosities radiated by 
the AGNs of the main progenitors in the upper panels of Fig. [T3] 
for our three simulated galaxies. A main progenitor is defined by 
applying the friends-of-friends algorithm for stellar particles with a 
linking length ln n ^ = 0.02 (I), where (I) is the the mean separation 
of dark matter particles. To estimate the luminosity of an AGN, we 
assume that the efficiency of conversion of rest mass energy of the 
accretion flow into radiation (Lagn/(Mc 2 )) is r/su = 10% for a 
standard thin disk and for a RIAF we use the equation 

r?RiAP = 0.1 ( -^- + <5 + KT 3 ) (26) 

where 8 repre sents the fractio n of the turbulent energy which heats 
the electrons dOuataert]|200ll) . While the value of 8 is quite uncer- 
tain, we e mploy 8 = 0.1 as com monly assumed in spectral models 
of RIAFs (Nemmen et al. 2006). Note that the radiative efficiency 
for a thin disk can be a factor of a few higher than 0.1 in the case 
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Figure 13. Upper panels: Bolometric luminosities radiated by the stars in the main progenitors (solid lines) and by the AGNs of the main progenitors (dotted 
lines). A 10% radiative efficiency is assumed for standard thin disks and the efficiency given by eq. {26} is assumed for RIAFs. Lower panel: Raw feedback 
luminosities ejected by SNe (II and la) in the main progenitors (solid lines) and by jets from the central BHs (AGNs) of the main progenitors (dotted lines). 
From left to right, the 'no-fb', 'weak-fb', and 'reference' simulations are shown. 



of a spinning BH, and may reach efficiencie s as high as 42% for an 
extreme Kerr BH jNovikov & Th orne 1973). 

It is found that the main progenitors are always brighter than 
the AGNs for these galaxies as we expect for disk-dominated galax- 
ies. Since the radiation drag is more efficient for more compact star 
forming region, we expect that the significance of the AGN lumi- 
nosity would increase for galaxies which are undergoing gas-rich 
mergers. There is a notable difference among our models when the 
accretion flows become RIAFs after z ~ 1. While all the AGNs 
have much lower luminosities than their host galaxies at low red- 
shift because of the low accretion rates and the low radiative ef- 
ficiency of the RIAF, only the AGN in the 'reference' galaxy has 
almost zero luminosity right after the accretion disc switches to a 
RIAF. In this simulation, the low accretion rate, due to the suppres- 
sion of star formation around the BH by AGN feedback, makes the 



AGN invisible most of the time. In the 'no-fb' and 'weak-fb' simu- 
lations, the AGNs still continuously emit light for a few giga years 
after the onset of the RIAFs. This result suggests that if 10% of the 
jet luminosity is used to heat the ambient halo gas AGNs become 
quite dim and can be hardly observed after the starburst regime. It 
should be noted that it is beyond the scope of this work to incorpo- 
rate in the simulations possible inclination effects in the appearance 
of the AGN, caused by dust obscuration inside the nucleus as in the 
obscuring dust torus model i nvoked in unified models of AGN (e.g. 
I Antonuccill 1 9931: iKrolikl 1 999h . 

In the lower panels of Fig.[T3]we show the raw feedback lu- 
minosities ejected by the total SNe in the main progenitors and jet 
luminosities produced by the AGNs computed by using equations 
( 1201 ) and ( 12 1 1 . Until the end of the starburst, the feedback energy is 
dominated by energy produced by SNe except for a few temporary 
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RIAF modes seen in the 'weak-fb' simulation. The jet luminosi- 
ties suddenly exceed the SN luminosities when the accretion flow 
becomes radiatively inefficient. Since in the 'no-fb' and 'weak-fb' 
simulations, the efficiencies, jjagn by which the AGN jet energy 
is given to the surrounding halo gas are small (?7agn = and 
0.01 respectively), the effects of the AGN feedback are none or 
negligible in these simulations. In the 'reference' simulation, the 
efficiency, 7/agn = 0.1 is sufficiently high to suppress the star for- 
mation around the BH therefore quenching the mass accretion to 
the BH. This is an important feature of our AGN feedback. The im- 
pulsive events in Fig.ll3lrepresent epochs where the AGN would be 
seen as radio-loud, and they constitute the radio-loud mode of the 
AGN. The duration in which the AGN is in the radio-loud mode 
is strongly dependent on the feedback efficiency, since the AGN 
feedback terminates the radio-loud mode through the suppression 
of gas cooling and hence star formation. Interestingly enough, in 
the 'reference' simulation, the AGN spends most of the time in the 
radio-quiet mode and only a fraction 8.8% of the time in the radio- 
loud mode. Contrarily, in the 'no-fb' and 'weak-fb' galaxies, the 
AGNs spend ~ 25% of the time in the radio-loud mode. Note that 
we are likely to overestimate the star formation rate at the centre 
and hence AGN activity particularly in the quiescent star forming 
regime because of the angular momentum transfer due to the artifi- 
cial an d numerical viscosities which can bring gas particles to the 
centre l lOkamoto et alj|2003l ; lkaufmann et"ail 20071) . 

5.3 Dependence on feedback parameters and 
implementations of feedback 

In §5.1, we compared simulations with different AGN feedback ef- 
ficiencies. In this sections, we compare simulations which all have 
strong AGN feedback, but which differ in other aspects. We com- 
pare three additional models. The models investigate the coupling 
of the AGN feedback, the role of BH spin and the effect of the wind 
speed generated by SN feedback. 

In the first model, AGN feedback energy is not given to the 
ambient halo gas (diffuse gas) but deposited to the nearest 40 ISM 
particles (dense gas). We use 77AGN = 1 so that all the jet en- 
ergy is used to heat the surrounding ISM because a model with 
tjAGN =0.1 does not show any visible difference from the 'no-fb' 
model. We call this model the 'ism-fb' model. In the second model, 
we assume the higher spin, j — 0.9, and the maximum feedback 
efficiency, tjagn = 1, as we may need such strong feedback in 
order to explain the X-ray properties of clusters of galaxies (e.g. 
Nemmen et al. 2007). As a result, the efficiency of AGN feedback 
in this model becomes 17 times higher than in the 'reference' model 
. We dub this model the 'max-fb' model. In the third model, we use 
the same AGN feedback as in the 'reference' model but employ a 
slow wind speed, 250 km s _1 for SN feedback in order to explore 
effects of the wind speed. We refer this model as the 'slow-wind' 
model. 

In Fig. [14] we show the edge-on and face-on views of the stel- 
lar distributions of the galaxies produced by the three additional 
models. The 'ism-fb' model is similar to the 'reference' model. 
As expected the disk component of the 'max-fb' galaxy is less 
significant than that of the 'reference' model since stronger AGN 
feedback suppresses the disk formation more strongly. The 'slow- 
wind' model exhibits completely different morphology from oth- 
ers. There is no sign of a disk component. 

Now we carry out quantitative comparison by showing the star 
formation histories of the galaxies and the accretion histories to the 
central BHs (Fig.ll5t. In spite of the fact that we employ an order of 



Figure 14. Same as Fig. \5\ but for the 'ism-fb' (left), 'max-fb' (middle), 
and 'slow-wind' (left) models. 

magnitude higher AGN feedback efficiency in 'ism-fb' model, the 
star formation history of this galaxy is very similar to that of the 
'reference' model. It suggests that delivering AGN feedback en- 
ergy to the diffuse hot gas causes the stronger feedback effect than 
depositing it into the surrounding dense ISM where cooling time is 
very short. The 'max-fb' galaxy shows a resembling star formation 
history to the 'reference' model (upper-middle panel) at high red- 
shift (t < 6 Gyr) where the AGN harbours a thin accretion disc. 
Once the accretion flow becomes a RIAF, strong AGN feedback 
suppresses the star formation dramatically and this galaxy have far 
less star formation after t ~ 6 Gyr than the 'reference' galaxy. The 
'slow- wind' galaxy has a completely different star formation his- 
tory from others. Owing to the slow wind speed, it has a higher star 
formation rate at high redshift including two burst events at t ~ 2 
and 3 Gyr. The AGN in this galaxy becomes underfed earlier than 
in other models (at t ~ 4 Gyr) and the AGN feedback strongly 
suppresses the subsequent star formation. We also show the model 
which employs the same wind speed as the 'slow-wind' model but 
without AGN feedback (dashed line in the upper-right panel). By 
comparing these two slow-wind models, we can see that the AGN 
feedback is responsible for the suppression of the star formation 
at low redshift. The bulge-to-total mass ratio of the 'slow-wind' 
model is 0.9 and the u — r colour of this galaxy is consistent with 
the mean colour of the red population at z = 0. 

The evolution of the accretion rates in the 'ism-fb' and 'max- 
fb' models are similar to that in the 'reference' model particu- 
larly at high redshift (t < 6 Gyr) where the accretion disc is thin 
(m > 10~ 2 ). In these two models, the accretion flows become RI- 
AFs after the starbursts as in the 'reference' model. After that the 
mass accretion onto the central BH in the 'ism-fb' model is more 
strongly suppressed than in the 'reference' model, while the 'ism- 
fb' model has a higher star formation rate. Since the AGN feedback 
energy is deposited into the surrounding star forming ISM particles 
in this model, it can terminate star formation around the BH more 
directly than in the 'reference' model where the AGN feedback en- 
ergy is distributed to the halo gas. In the max-fb model, there is no 
accretion event after t ~ 6 Gyr. The strong AGN feedback in this 
model removes all the low angular momentum gas around the BH 
and prohibits subsequent star formation around the galactic centre. 
The 'slow-wind' model shows a different accretion history from 
other models. Since the starbursts occur earlier in this model ow- 
ing to the slow wind speed, the AGN also becomes underfed earlier 
(t ~ 4 Gyr) than in other models. Contrary to other models, the 
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Figure 15. Upper panel: star formation histories for the 'ism-fb' (left), 'max-fb' (middle), and 'slow-wind' (right) models. The star formation history for the 
reference model (dashed line in Fig. ref'SFH) is indicated by the dotted line. In the right panel, we also show the 'slow-nofb' model in which we use the same 
wind speed as in the 'slow-wind' model and do not include AGN feedback (dashed line) for comparison. Lower panel: same as the lower panel of Fig.[8]but 
for the 'ism-fb' (left), 'max-fb' (middle), and 'slow-wind' (right) models. 



mass accretion continues at t ^ 4 Gyr. This is because the slow 
wind speed, 250 km s _1 , does not allow gas particles to escape 
from galactic centre and the AGN feedback energy is not given 
to the star forming gas. The accretion rate might become more 
episodic if we deposit a significant fraction of the AGN feedback 
energy to the central star forming gas as in the 'ism-fb' model. 

In Fig.|T6] we show the relations between mass of the central 
BHs of the main progenitors and the mass of the host spheroids at 
z = 3, 2, 1.5, 1, 0.5, 0.2 and 0. The relations for 'ism-fb' and 'max- 
fb' models evolve similarly to the relations for the models pre- 
sented in Fig. [TJJ The relation for the 'slow-wind' model reaches 
the local relation earlier (at z = 2) because this galaxy undergoes 
starbursts earlier than other galaxies. At low redshift, this model 
has slightly smaller BH mass to spheroid mass ratio than in other 
models. This is because infalling satellites increase the mass of the 
spheroid without changing the mass of the central BH very much 
when these infalling galaxies are disk dominated. Due to the slow 
wind speed, these infalling systems have more stars than in other 
models. Fig. [TT] and Fig. [16] tell us an important feature of the ra- 
diation drag model, that is, the Mbh — Msp relation after the last 



starburst does not depend on either the wind speed or efficiency of 
AGN feedback. In all models, the central BH masses are ~ 0.1% 



of the mass of their host spheroids, Mbh 
accretion becomes radiatively inefficient. 



10 M S p, when the 



COMPARISON WITH OBSERVATIONS OF AGN IN 
SPIRAL GALAXIES 



As we have shown in S|5.2| for z < 1 the bolometric luminos- 
ity of the AGN is feeble compared to the host galaxy and the nu- 
cleus undergoes several recurring episodes of activity. During these 
episodes when the nucleus is active, its characteristic accretion rate 
is in the range m w 10 -4 — 10 -3 , while the typical nuclear bolo- 
metric luminosities and jet powers are Lboi ~ 10 40 — 10 42 



and L v 



erg s 



jjet ~ 10 —10 erg s~ respectively. Given these prop- 
erties, at z < 1 during which the galactic disk is formed, the 
AGN can be considered to be in a highly sub-Eddington, low- 
luminosity or Seyfert state which undergoes a stochastic cycle of 
activity. The Seyfert nucleus is sometimes turned on, but most of 
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Figure 16. Same as Fig. 1 1 1 1 but for the 'ism-fb' (pluses), 'max-fb' (dia- 
monds), and 'slow-wind' (triangles) models. 



the time is "off" or inactive with m ~ 0, resembling the nucleus 
of a normal spiral galaxy. The "on" episodes are induced by the 
inward motion of clumps of cold gas into the centre of the galaxy, 
which provide a replenishment of new stars. The radiation drag or 
Poynting-Robertson effect caused by stellar radiation in the cen- 
tral star-forming region then dissipates the angular momentum of 
the gas and fuels the AGN. Our model suggests that reasonably 
powerful jet outbursts accompany each "on" episode of the low- 
luminosity AGN, suggesting that, in analogy to the observed states 
of AGNs, during the active phases the Seyfert nucleus would be 
seen as radio-loud, becoming radio-quiet during the "off" phase. In 
this section, we compare the properties of our simulated AGN with 
those of observed ones in disk galaxies. 

In the reference model, the range of values obtained for the 
AGN bolometric luminosity for z < 1 agrees well with the median 
observational values for LINERs (median Lboi ~ 10 41 erg s _1 ) 
and Seyferts (median Lboi ~ 10 42 erg s _1 ) i nferred fr om the 
Palomar spectroscopic survey of nearby galaxies (Ho 2004). Since 
the Seyfert episodes of the central AGN coincide with jet out- 
burst events during which the AGN would be seen as radio-loud 
in our simulations, this suggests that observed Seyfert/LINER nu- 
clei would be mostly radio-loud objects rather than radio-quiet 
given their formation history, contrar y to the gene r ally held no - 
tion existent in the literat ure (e.g., Krolik 19991: Laorl 2000 ) . 
iHo & Ulvestadl l l200lh and iTerashima & Wilsorj ( 120617 however 
investigated several low-luminosity AGN and their radio-loudness 
classification, obtaining that a large fraction of Seyfert/LINER nu- 
clei ar e radio-loud, a finding that corroborates our results. In par- 
ticular, |Ho|& Ulvestad feOOll) obtained that > 60% of the sources 
in their sample are radio-loud. 

In our simulations, the numerical resolution is not sufficient to 
resolve the details of the propagation of the jet as it flows out of the 
AGN and how it transfer its power to the surrounding medium, for 
instance the spatial resolution is 0.5 kpc. As such, we assume that a 
fraction of the jet power is injected as thermal energy in the diffuse 
gas around the AGN on kpc-scales. The size of the region which 



is affected by the feedback from the AGN may be as large as w 
10 — 15 kpc, although the jet itself is not resolved in the simulation. 
Two important questions arise: (a) are the heating effects due to 
jets generated in Seyfert galaxies able to affect kpc-scales and (b) 
can they be a source of feedback? Observations provides appealing 
clues to answering these questions. 

Several observational studies of Seyferts suggest that Seyfert 
jets may indeed reach kpc-sc a les an d be yond. The Very Large 
Array surveys of IColbert etall d 19961) and lballimore et al'j (2006) 
find that extended kpc-scale radio outflows most probably asso- 
ciated with jets are common in Seyfert galaxies, although they 
cannot rule out entirely the possibility that these outflows arise 
from starburst-driven winds in some sources. These radio outflows 
have random orientations with respect to the host galaxy, a result 
which is maintained for jets observed fr om sub-kpc to kpc scales 
dKinnev et al]l2000l : ISchmitt et alj|200ll) . and have distorted, un- 
collimated morphologies. Therefore the observed radio outflows 
could provide a broad effective surface for depositing the jet power 
in the surrounding gas, as envisioned in our simulations. Several 
targeted observations of individual Seyfert galaxies also demon- 
strate thatkjjc-scjilej^id^ ie^^ 



1998; [Cecil eTaL 


2000; Whittle & Wilsorj|2004l; iKeel et al.ll2006: 


Kharb et alj|2006 


.iMiddelbere et all2007l). 



X-ray observations also reveal kpc-sc ale outflows of X- 
ray-emitting hot plasma in Seyfe rt galaxies dColbert et alj|l998l : 
iRupke. Veilleux. & Sanders! 120051) . cospatial with the large scale 
radio outflows. Two favoured interpretations are that these kpc- 
scale X-ray outflows are generated either by the AGN jets that en- 
train and heat gas as they flow out of the nucleus, or through heat- 
ing of the surrounding gas via thermalisation of the kinetic power 
of the jet, generating a wind. Either way, these processes bear re- 
semblance to the generation of cavities or bubbles of X-ray emit- 
ting plasma in the central elliptical gala xies of clusters of galaxies , 
which are inflated by the radio jets (e.g . , iDalla Vecchia et alj2004l ; 
lBirzanetal1l2004l ;l Fabia n et al. 2006). These observations there- 
fore provide evidence for the possible importance of feedback due 
to the radio jets in Seyfert galaxies, which could heat the circum- 
nuclear plasma and generate X-ray emitting outflows. 

Observational estimates of the kinetic power carried 
by the jets in Seyferts using different methods indicate 



that it is 



in 
in the 



range 



10 4 



erg_ 



(e.g. , 



Wilson & Ulvestad l983l;lMorganti Oosterloo. & Tsvetanovl 19981: 



iBicknell et ai]|l998l : ICapetti et al.lfl9991 ; iKharb et alj|2006t) . This 
observational range is in agreement with the values of the jet power 
which are generated during the radio-loud episodes for z < 1 in the 
simulation, L jct » 10 42 - 10 43 erg s" 1 (Fig.[[3]l. It is worth not- 
ing that the jet power in our simulation is estimated using a physical 
model for the energetics of the jet, which depends on fundamental 
parameters of central engine, especially the accretion rate which is 
self-consistently calculated. 

The fraction of the jet power generated by the AGN which is 
actually available to heat the gas on kpc-scales is a matter of debate. 
In the 'reference' simulation we adopt that a fraction jjagn = 0.1 
of the jet power is deposited as thermal energy of the diffuse gas, 
although it is possible that this may be an overestimate of the value 
of tjagn- For instance, based on ram pressure and terminal ve- 
locity co nsiderations of the lob es of kpc-scale radio outflows in 
Seyferts, Gallim ore et alj ((2006) argues that the Seyfert jets lose 
almost all their power within the inner kpc, implying that 77agn 
may be considerably smaller than 0.1. The other simulations we 
present encompass the uncertainty in the value of ^agn- In this re- 
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spect, the simulation adopting 77AGN = 0.01 does not show any 
significant difference from the simulation without AGN feedback 
(??agn = 0). On the other hand, when 100% of the jet energy is 
transferred to the surrounding star forming gas ('ism-fb' model), 
most of the feedback energy is rapidly radiated away, and the re- 
sults are very similar to the 'reference' model in which 10% of jet 
energy is given to the ambient halo gas. Thus, a better knowledge 
on the "microscopic" physics of jet-ISM interactions is required in 
order to improve the understanding of the effects of AGN feedback 
in Seyfert galaxies. 

Our simulations show that during the cosmological forma- 
tion of a typical spiral galaxy, a low-luminosity AGN (Seyfert 
or LINER) is formed which has a stochastic cycle of activity, 
with many short-lived recurrent accretion events accomp anied by 
jet ou tbursts. This finding agrees with the conclusion of Sanders 
dl984l) that the Seyfert activity is the result of short-lived stochas- 
tic accretion events. He arrived at this conclusion by consider- 
ing plausible mechanisms that may fuel the central black holes 
of Seyferts, which are stochastic in nature, and our simulations 
track the evolution of the AGN by implementing one of these 
mechanisms: the accretion of molecular clouds to the galactic cen- 
tre. Our findings therefore provide support to the assumption of 
Hopkins & Hernquistl J2006l) . central to their work, that the central 
supermassive BHs of Seyferts stochastically increase their mass by 
accretion of cold gas. 

In the period z < 1 after the galactic disk is formed, the AGN 
is active only during a fraction ~ 9% of the its life in the 'reference' 
model, « 25% in the 'weak-fb' model and ~ 1% in the 'ism- 
fb' model. This implies that the lifetime of Seyfert activity in the 
spiral galaxy is « 10 9 yr for the 'reference' or 'weak-fb' model 
and ~ 10 s yr for the 'ism-fb' model. This is in agreement with 
statistical Seyfert lifetime ~ 10 8 — 10 9 yr estimated by cou nting the 
fraction of spiral galaxies that are observed to be Se yferts dWoltierl 
ll968tlSandedl984lHo et all 19971) . I Sanders! ( 1 19841) estimated that 
each single episode of activity would not last longer than 10 6 yr 
based on the size of Seyfert jets, a prediction w hich is in agreement 
with the detailed study of iKharb etall fcOOd) . Unfortunately the 
time-step we used to output snapshot files is 16.8 Myr and therefore 
we may miss many of such short-lived events. Considering this, 
the 'ism-fb' model may be more comfortable than the 'reference' 
model because there is no long-lived accretion event in the 'ism- 
fb' model for z < 1 (lower-left panel of Fig.[T5l>. It is also possible 
that the observed short duration of activity is not caused by the 
large-scale mass accretion to the galactic centres but by the physical 
processes in the accretion disk (e.g., instabilities) which we do not 
take into account in this paper. 



7 SUMMARY AND CONCLUSIONS 

In this paper, we have introduced a new methodology that enables 
us to simultaneously model star formation, SN feedback, BH ac- 
cretion, and AGN feedback in cosmological simulations of galaxy 
formation. Our treatment of black hole accretion is based on a the- 
oretical model of the mass accretion into the galactic centre due to 
radiation drag. This connects the kpc-scale star formation activity 
in the galactic centre with the BH accretion ra te. Furthermo r e, mo- 
tiva ted by the recent se mi-analytic models of ICroton et al. ( 2006) 
and lBower et alj d2006l) . we have distinguished two distinct accre- 
tion modes of BH fuelling, namely standard (optically thick, geo- 
metrically thin) Shakura-Sunyaev disks and (geometrically thick, 
optically thin) RIAFs. This allows us to distinguish radio loud and 



radio quiet AGNs depending on their mass accretion rate. In this 
model, we only consider AGN radio feedback in which we assume 
a fraction of the jet power is thermally coupled with the surround- 
ing diffuse gas, and assume RIAFs are responsible for production 
of powerful jets. As the first application of this model, we have per- 
formed cosmological simulations of the formation of a disk galaxy 
and its coevolution with the central AGN. 

Our results reveal a fundamental AGN-starburst connection 
during the evolution of the galaxy and the central AGN. During the 
starburst period (z > 1) there is a lot of gas available to be accreted 
and the AGN is fuelled at high accretion rates, as a consequence the 
accretion disk is in the standard/thin state and the AGN is relatively 
luminous and produces weak jets (radio-quiet state). The starburst 
phase ceases at z ~ 1 because most of the gas at the galactic cen- 
tre is consumed by star formation or blown out by the associated 
SN feedback. Afterwards, the AGN fuelling rate drops and the ac- 
cretion disk then switches to a radiatively inefficient (RIAF) state. 
Due to a combination of small accretion rates and the low radiative 
efficiency of the RIAF, the AGN becomes a low-luminosity AGN 
or Seyfert nucleus which is almost invisible compared to the host 
galaxy. At this point even a little star formation activity around the 
central BH triggers intense jet production that regulates further gas 
cooling to the centre. Thus, in the quiescent star forming regime 
after z ~ 1 the AGN has a stochastic fuelling, undergoing several 
short-lived recurrent episodes of activity. During these accretion 
episodes the AGN brightness increases and jet outbursts are gen- 
erated, which correspond to the radio-loud state. In-between these 
Seyfert episodes the nucleus is turned off and the galaxy is inactive. 

Below z ~ 1 AGN feedback becomes important if 10% of jet 
energy is available to heat the ambient diffuse halo gas (or 100% is 
fed into the surrounding ISM). This powerful episodic AGN feed- 
back suppresses star formation in the disc and almost completely 
halts the star formation in the bulge. This effect is minor when the 
wind speed due to SN feedback is fast compared to the circular 
velocity of the host halo because the SN feedback alone can sup- 
press the star formation in the bulge reasonably well. When the 
wind speed is slow compared to the circular velocity of the host 
halo, the AGN feedback becomes more important. This leads to the 
formation of an old, red elliptical galaxy in the 'slow-wind' model. 

Many properties of the simulated low-luminosity AGN and 
its jets are in broad agreement with observations of Seyfert galax- 
ies and their jets: namely its typical nuclear bolometric luminosity 
and radio luminosity, the coincidence of the radio-loud state with 
the Seyfert episodes, the fact that the jet may reach kpc-scales and 
strongly affect the surrounding gas on these scales, and the pre- 
dicted AGN lifetime and duration of individual episodes of activ- 
ity. The stochastic pattern of nuclear activity th at emerge s from our 
simulation seems to confirm the prediction of Sanders] d 19841) for 
Seyfert galaxies. 

The ratios between the central BH mass and the mass of the 
host spheroid at 2 = in our simulations are ~ 10 -3 regard- 
less of the strength of either SN or AGN feedback. This result is 
direct outcome of the radiation drag model which relates the star 
formation rate in the central star-forming region and the mass ac- 
cretion rate to the central BH. The ratio is in good agreement with 
the observed ratio dKormendv & Richstondl 19951; iMagorrian et al] 
1 19981 : iMerritt & Ferraresei l200ll ; iMcLure & Dunlod |2002|) , al- 
though the efficiency of the radiation drag we assumed to ob- 
tain this ratio is somewhat higher than the value suggested by 
Kawakatu &UmemuriJj2002l) . However, as we have discussed, the 
geometry and density structure of star forming regions can easily 
change the efficiency. We also note that recent semi-analytic mod- 
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els ( iBaugh eta! 20051; Nagashima et a l. 2005a. b) and cosmologi- 
cal simulations (OEFJ05) suggest that stars may be born with top- 
heavy IMFs (and thus stronger radiation fields and higher radiation 
drag efficiencies) in starbursts. The accretion rate by radiation drag 
can be much higher with this assumption because stellar popula- 
tions with top-heavy IMFs produce more radiation. 

The self-consistent treatment of BH accretion and associated 
AGN feedback in galaxy formation enables us to explore the in- 
terplay between galaxies and AGNs in cosmological simulations 
of galaxy formation. The self-regulation of mass accretion due to 
feedback from RIAFs produces galaxies whose AGN properties are 
well matched to the observations. The model is also qualitatively 
consistent with recent semi-analytic models that successfully re- 
produce many of the global properties of the galaxy population. 
The effects of our AGN feedback implementation on formation 
of groups of galaxies, where more powerful AGN feedback from 
larger BHs is expected, will be presented in a forthcoming paper. 
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